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Abstract 

A program to generate codes in Fortran and C of the full Magne¬ 
tohydrodynamic equations is shown. The program used the free com¬ 
puter algebra system software REDUCE. This software has a package 
called EXCALC, which is an exterior calculus program. The advan¬ 
tage of this program is that it can be modified to include another 
complex metric or spacetime. The output of this program is modified 
by means of a LINUX script which creates a new REDUCE program 
to manipulate the MHD equations to obtain a code that can be used 
as a seed for a MHD code for numerical applications. As an example, 
we present part of output of our programs for Cartesian coordinates 
and how to do the discretization. 


1 Introduction 

Nowadays, there are programs that generate codes for a given mesh 
arrangement, but there is no a program to generate the differential 
equation in a given space-time. Now, we show how it is possible by 
means of a EXCALC package [16] of REDUCE [8] and a LINUX script. 
REDUCE is a free computer algebra system (CAS) software intended 
to algebraic manipulation. In that sense, it is similar to Mathematica 
or Maple. The EXCALC package solves problems using the Cartan 
formalism or differential forms [1, 7]. 
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Due to the advent of new technology, there is an interest in plasma 
simulations and visualizations e.g. [4, 6, 14], There are a lot of infor¬ 
mation on these subjects in the Internet [2, 12, 9, 10, 17]. Moreover, 
the mesh or grid generation is important in designing complex engi¬ 
neering structures or complex nonsymmetric systems [15]. 

We wrote a REDUCE program (MHD.red) and a long UNIX script 
(discretized-mhd) mainly for Magnetohydrodynamics (MHD), but 
the program and script we present here are easy to modify for another 
purposes. Our program could be used together with a mesh generation 
code to solve complex problems. 

The paper is organized as follows. In section 2, a brief description of 
MHD is presented. A short summary of EXCALC is given in section 
3. The explanation of the program is described in section 4. Moreover, 
we present the results of the discretization to Cartesian coordinates as 
an example. The conclusions and future work is discussed in section 
5. 

2 MHD 

MHD is the study of the dynamics of electrically conducting fluids. 
Examples of such fluids include plasmas, liquid metals, and salt wa¬ 
ter or electrolytes. Plasmas can be regarded as fluids since the mean 
free paths for collisions between the electrons and ions are macroscop- 
ically long. Thus, collective interactions between large numbers of 
plasma particles can isotropize the particles velocity distributions in 
some local mean reference frame, thereby making it sensible to de¬ 
scribe the plasma macroscopically by a mean density, velocity, and 
pressure. These mean quantities can then be shown to obey the same 
conservation laws of mass, momentum and energy for fluids. 

The fundamental concept behind MHD is that magnetic fields can 
induce currents in a moving conductive fluid, which in turn creates 
forces on the fluid and also changes the magnetic field itself. The set 
of equations which describe MHD are a combination of the Navier- 
Stokes equations of fluid dynamics and Maxwell’s equations of elec¬ 
tromagnetism. 

An important physical parameter that defines the way MHD treats 
the systems is resistivity. When the resistivity is very low plasmas can 
be treated as perfect conductors (ideal MHD in the infinite magnetic 
Reynolds number limit) and the the lines of force appear to be dragged 
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along with the conductor. 

When the resistivity cannot be neglected the magnetic field can gener¬ 
ally move through the fluid following a diffusion law with the resistivity 
of the plasma serving as a diffusion constant. According to this the 
solutions to the ideal MHD equations are only applicable for a limited 
time before diffusion becomes too important to ignore. The diffusion 
time across a solar active region is hundreds to thousands of years. 
The interested reader may consult the references [3, 5, 11, 13]. 

There are many applications of MHD but in this article we are inter¬ 
ested in those related to astrophysics and cosmology. Plasma is the 
constituent of many astrophysical objects of the Universe: stars, neb¬ 
ulae, the interplanetary, the interstellar and the intergalactic media. 
Many astrophysical problems require the treatment of magnetized flu¬ 
ids in dynamical, strongly curved spacetimes. Such problems include 
the origin of gamma-ray bursts, magnetic braking of differential ro¬ 
tation in incipient neutron stars arising from stellar core collapse or 
binary neutron star merger, the formation of jets and magnetized disks 
around newborn black holes and many others [18]. 

MHD Equations: 

Continuity Equation 

The conservation of mass in a fluid tell us that 



( 1 ) 


where p is the density, v is the fluid velocity and = pv. 


Maxwell Equations 

The electric Gaufi law is 


(2) 


eoV • E = pe 


where E is the electric field, pe is the electric charge density and cq is 
the vacuum permittivity. 

The magnetic Gaufi law is 


V • B = 0. 


(3) 


where B is the magnetic field. 
The Ampere law is 
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V X B = fioJ + fioeo — . (4) 

where J is the electric current density and /io is the vacuum perme¬ 
ability. 

The Faraday law is 


V X E = - 


SB 


(5) 


Equation of Motion 

The equation of motion of a fluid or the Navier-Stokes equation with 
electromagnetic fields is 

(9T 

y = Fi + V.[P-pvv]. (6) 

where P is a the stress tensor with components 

Pij = ~P^ij + Pij (7) 

P is the pressure and Tij represents the viscosity stress tensor compo¬ 
nents, which are given by 

_ (dvi dvj 2 \ „ e 

+ ^ - 3^'® 

with ^ and v are the first and second viscosity coeficients. 

Fl is the volume Lorentz force given by 


Fi = peE -I- J X B 
= V • T — eo/^o 


(9) 


'di 


where S is the Poynting vector 
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S = — E X B. 

Ato 

and T is the electromagnetic tensor with following components 


Tij — ^QpiPj H- BiBj — I Co 

/^o 


, B^\ 
2/^0 / 








( 10 ) 


( 11 ) 
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The Ohm Law 

The Ohm law is 


E = ryj — V X B, 

ry is the resistivity. 

Equation of State 

The equation of state for an ideal gas is 

P = pRT. 

where T is the temperature and R is the gas constant. 


Equation of Energy Conservation 

The equation of energy conservation is 


where 


m 



[V • Q - J • E]. 


Q 



Js + q ~ P 


V 


and e is the internal energy 


e = CvT 

here 7 is the adiabatic exponent 


P 

( 7 - l)/5’ 



and q is the heat flux vector 


R + Cv 
Cv 


q = —kVT, 

with k is the thermal conductivity. 


( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 
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3 EXCALC 


EXCALC is a package under REDUCE intended to solve problems 
with Cartan exterior calculus. The advantage of using it, is that given 
equations can be written independent on the coordinate system. The 
outcomes of the calculations from these equations are expressed in the 
chosen coordinate system. Here, it is a list of the most used commands; 

Exterior multiplication 
Partial differentiation 
Hodge * operator 
Inner prodnct 
Lie derivative 
Declaration of a coframe 
Exterior differentiation 
Declaration of implicit dependencies 
Declares the frame dual to the coframe 
Clause of CDFRAME to specify a metric 
Declaration of exterior forms 
The principal operations with 1-forms or vectors are listed below 

Cross product ^(v^B) (B and v are vectors or 1-forms) 

Nabla d P (P is a function or scalar) 

Divergence #d^ F (F is a vector or 1-form) 

Curl F (F is a vector or 1-form) 

Eor more information about this REDUCE package, the interested 
user may consult the references. 


0 

# 


COFRAME 

d 

FDOMAIN 

FRAME 

METRIC 

PFDRM 


4 The Program 

The REDUCE program MHD.red to transform the MHD equations 
into Cartesian coordinates is included as an appendix. The script 
discretized-mhd and the generating code are not include in this pa¬ 
per, but it can be emailed send by the authors if requested, and in 
a near future will be available at our Space Research Center Web¬ 
page http://cinespa.ucr.ac.cr/software/. 
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4.1 Example of the Script for Cartesian Coor¬ 
dinates 

Now, we will present the modified output of the programs for a given 
MHD equations using an UNIX script. Let us consider the induction 
equation: 


where E = 77 — v x B and //oJm = V x B, without taking into 
account the time variation of the electric field. 

For example, the outcome from our program for the y component of 
the induction equation (Cartesian coordinates) is 
induct( - x2) := @(b( - x2),t) 

- b( - xl)*@(v( - x2),xl) + b( - x2)*@(v( - xl),xl) 

+ b( - x2)*@(v( - x3),x3) - b( - x3)*@(v( - x2),x3) 

- @(b( - xl),xl)*v( - x2) + @(b( - x2),xl)*v( - xl) 

+ @(b( - x2),x3)*v( - x3) - @(b( - x3),x3)*v( - x2) 

+ (@(b( - xl),xl,x2)*eta - @(b( - x2),xl,xl)*eta 

+ @(b( - xl),x2)*@(eta,xl) - @(b( - x2),xl)*@(eta,xl) 

- @(b( - x2),x3)*@(eta,x3) + @(b( - x3),x2)*@(eta,x3) 

- @(b( - x2),x3,x3)*eta + @(b( - x3),x2,x3)*eta)/mu_0 


As this output is not easy to handle, we translate such equation to a 
FORTRAN or C code, then we wrote a script that modify this outcome 
into equations with a prescribed discretization method. To this aim, 
we employ the UNIX command sed [19]. This command change the 
output using the prescribed discretization method. For instance, the 
following commands 

sed -e ’s:@(b( - xl),t):(bx(n+l,i,j,k)-bx(n,i,jjk))/dt:g’ 
mhdeqs r mhdeqsnew 

sed -e ’s:@(v( - x 2 ),xl):(vy(n,i+ 1 ,j,k)-vy(n,i- 1 ,jjk) 
)/( 2 *dx):g’ mhdeqs r mhdeqsnew 


will produce the translation 

dhx b^{n + 1 , i, j, k) - b^{n, i, j, k) 
dt At 

dvy Vxjn, i + 1 , i, k) - Vxjn, i - 1 , j, k) 
dx ^ 2Ax 
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The output of our REDUCE program is mhdeqs (see Appendix) and 
the mhdeqsnew, which should be created before one runs the script, 
will contain all modifications of the entire system of equations. 

4.2 Example of the Code Generation 

Here, we show a part of the output after the script is run. In this case, 
the y component of the induction equation is converted into 
induct_y := (by(n+1,i,j,k)-by(n,i,j,k))/dt 
+ ((bx(n,i+1,j+1,k)-bx(n,i+1,j-1,k) 

-bx(n,i-1,j+1,k)+bx(n,i-1,j-1,k))/(4*dx*dy)*eta(n,i,j,k) 

+ (bx(n,i,j+l,k)-bx(n,i,j-l,k))/(2*dy)*(eta(n,i+l,j,k) 
-eta(n,i-l,j,k))/(2*dx) 

- (by(n,i+2,j,k)-2*by(n,i,j,k)+by(n,i-2,j,k))/(4*dx*dx)* 
eta(n,i,j,k) 

- (by(n,i+l,j,k)-by(n,i-l,j,k))/(2*dx)*(eta(n,i+l,j,k) 
-eta(n,i-l,j,k))/(2*dx) 

- (by(n,i,j,k+2)-2*by(n,i,j,k)+by(n,i,j,k-2))/(4*dz*dz)* 
eta(n,i,j,k) 

- (by(n,i,j,k+l)-by(n,i,j,k-l))/(2*dz)*(eta(n,i,j,k+l) 
-eta(n,i,j,k-l))/(2*dz) 

+ (bz(n,i,j+l,k+l)-bz(n,i,j+l,k-l) 

-bz(n,i,j-1,k+l)+bz(n,i,j-l,k+l))/(4*dy*dz)*eta(n,i,j,k) 

+ (bz(n,i,j+l,k)-bz(n,i,j-l,k))/(2*dy)*(eta(n,i,j,k+l) 
-eta(n,i,j,k-l))/(2*dz))/mu_0 

Thus, it is possible to use this output as a FORTRAN or C code. 

5 Conclusion 

In this paper, we show how to produce quickly a computer code useful 
to implement simulations or visualizations. The automatic generation 
of the computer code is possible by means of REDUCE programs and 
UNIX scripts that modify and manipulate the output of the programs 
and generate computer codes with a given discretization method. The 
programs shown here, are flexible and are easy to modify for other 
purposes. 

Due to the rapid advance of technology, there is an interest in com¬ 
puter simulations and visualizations, in particular in plasma physics. 


Our program together with recent developed grid generation codes 
helps solve complex plasma phenomena. 


6 Appendix 

7o REDUCE Program MHD.red 

% It is based on an example of E. Schruefer 
out mhdeqs$ 
load excalc$ 
off nat$ 

% Problem: 

7 

7 Calculate in spherical coordinates system the MHD 
7 equations. 

7 coframe e r = d r, e theta = r*d theta, 

7 e phi = r*sin(theta)*d phi$ 

7 frame x$ 

7 fdomain v = v(t, r, theta, phi), 

7 p = p(t, r, theta, phi)$ 

7 

7 Calculate in cartesian coordinates system the MHD 
7 equations. 

coframe e xl = d xl, e x2 = d x2, e x3 = d x3$ 
frame x$ 

fdomain v=v(t,xl,x2,x3), p=p(t,xl,x2,x3)$ 
fdomain rho=rho(t,xl,x2,x3), je=je(t,xl,x2,x3)$ 
fdomain pstar=pstar(t,xl,x2,x3), b=b(t,xl,x2,x3)$ 
fdomain eta=eta(t,xl,x2,x3), energy=energy(t,xl,x2,x3)$ 
fdomain jm=jm(t,xl,x2,x3), jmp=jmp(t,xl,x2,x3)$ 
fdomain ee=ee(t,xl,x2,x3), a=a(t,xl,x2,x3)$ 
pform v(k)=0, je(k)=0, p=0, rho=0$ 
v := v(-k) * e(k); 

7 je := rho * v; 
je := je(-k) * e(k); 

7factor e, @$ 
factor e$ 
on rat$ 

7 

7 First we calculate the continuity equation. 
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7o drho/dt + nabla.je = 0 

pform conteq=0$ 

conteq := @(rho,t) + #d# je; 

7 

7 Next we calculate the equation of motion. 

7 d je/dt + nabla pstar + (nabla.je) v + je.nabla v 
7 + (nabla.b) b + b.nabla b = 0 
7 nabla.b = 0 

pform moveq=l, moveq(-k)=0, pstar=0, b(-k)=0$ 
b := b(-k) * e(k); 

7 pstar := p+(b(k) * b(-k))/(2*mu_0) ; 

7 moveq := @(je,t) + d pstar + (#d# je)*v 
7 + rho*((v(k) * x(-k)) |_ v - (l/2)*d (v(k) * v(-k))) 
7 + (#d# b)*b + ((b(k) * x(-k)) |_ b 

7 - (1/2)*d (b(k) * b(-k))); 

moveq := @(je,t) + d pstar + (#d# je)*v 
+ rho*((v(k) * x(-k)) |_ v - (l/2)*d (v(k) * v(-k))) 

+ ((b(k) * x(-k)) |_ b - (l/2)*d (b(k) * b(-k))); 
moveq(-k) := x(-k) _| moveq; 

7 

7 Next we calculate Gauss equation for the magnetic 
7 field. 

7 nabla.b = 0 
pform gauss=0; 
gauss := #d# b; 

7 

7 Next we calculate the induction equation, 
pform induct=l, induct(-k)=0, ee=l, ee(-k)=0, 
jm=l, jm(-k)=0, eta=0$ 

pform jmp=l, jmp(-k)=0, jmpp=l, jmpp(-k)=0, 
a=l, a(-k)=0, bb=l, bb(-k)=0$ 

pform induct2=l, induct2(-k)=0, ee2, ee2(-k)=0$ 
jm := (#d b)/mu_0; 
jm(-k) := x(-k) _| jm; 

7b := #d a; 
a := a(-k) * e(k); 
jmp:= jmp(-k) * e(k); 
bb := #d a; 

jmpp := (#d bb)/mu_0 - jmp; 
jmpp(-k) := x(-k) _| jmpp; 
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ee := eta=t=jm - ttCv^b); 
ee(-k) := x(-k) _| ee; 
ee2 := eta*jm - #(v~bb); 
ee2(-k) := x(-k) _| ee2; 
induct := @(b,t) + #d ee; 
induct(-k) := x(-k) _| induct; 

induct2 := @(a,t) + eta*(#d bb)/mu_0 - #(v~bb); 
induct2(-k) := x(-k) _| induct2; 

7o 

7 Finally we calculate the energy conservation. 

7 energy := rho * v**2/2 + p/(gamma-1) + b**2/(2*mu_0); 
pform energy=0, enconser=0$ 

7 energy := rho * (v(k) * v(-k))/2 + p/(gamma-1) 

7 + (b(k) * b(-k))/(2*mu_0) ; 

enconser := ©(energy,t) + #d# ((energy+pstar)*v 
- ((v(k) * x(-k)) _| b)*b); 

7 

clear v, x, b, p, pstar, eta, conteq, moveq, gauss, 
ee, energy, enconser$ 
remfac e, @$ 

remfdomain p, v, b, energy, eta, pstar, rho, je$ 

shut mhdeqs$ 

bye; 
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